Decomposition, Reformulation, and Diving 
in University Course Timetabling 



Edmund K. Burke, Jakub Marecek, Andrew J. Parkes 

q ! School of Computer Science, University of Nottingham 

O ' Jubilee Campus, Nottingham NG8 IBB, UK 

<n : 

c3 ■ Hana Rudova 

§ : 

^ ■ Masaryk University Faculty of Informatics 

Botanickd 68a, Brno 602 00, The Czech Republic 

XS1 

O ■ Abstract 
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In many real-life optimisation problems, there are multiple interacting components 
in a solution. For example, different components might specify assignments to dif- 
ferent kinds of resource. Often, each component is associated with different sets 
^ ■ of soft constraints, and so with different measures of soft constraint violation. The 

Q\ . goal is then to minimise a linear combination of such measures. This paper studies 

an approach to such problems, which can be thought of as multiphase exploitation 
of multiple objective-/value-restricted submodels. In this approach, only one com- 
putationally difficult component of a problem and the associated subset of objec- 
Q\ \ fives is considered at first. This produces partial solutions, which define interesting 

neighbourhoods in the search space of the complete problem. Often, it is possible 
to pick the initial component so that variable aggregation can be performed at the 
first stage, and the neighbourhoods to be explored next are guaranteed to contain 
feasible solutions. Using integer programming, it is then easy to implement heuris- 
tics producing solutions with bounds on their quality. 

Our study is performed on a university course timetabling problem used in 
the 2007 International Timetabling Competition, also known as the Udine Course 
Timetabling Problem. The goal is to find an assignment of events to periods and 
rooms, so that the assignment of events to periods is a feasible bounded colouring 
of an associated conflict graph and the linear combination of the numbers of viola- 
tions of four soft constraints is minimised. In the proposed heuristic, an objective- 
restricted neighbourhood generator produces assignments of periods to events, 
with decreasing numbers of violations of two period-related soft constraints. Those 
are relaxed into assignments of events to days, which define neighbourhoods that 
are easier to search with respect to all four soft constraints. Integer programming 
formulations for all subproblems are given and evaluated using ILOG CPLEX 11. 
The wider applicability of this approach is analysed and discussed. 
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1 Introduction 

1.1 Motivation 

There has recently been considerable progress in the development of mixed 
integer programming solvers (Nemhauser and Wolsey, 1988; Johnson et al., 
2000; Bixby et al., 2004). To a considerable extent, this progress is due to the 
introduction of new primal and improvement heuristics (Fischetti and Lodi, 
2003; Fischetti et al., 2005; Danna et al, 2005; Berthold, 2007a). These meth- 
ods can be said to "dive" into a particular region of the search space in order 
to explore it intensely. This is intended to restrict and simplify the search 
space so that the sub-instance is easier to solve. Currently, diving heuris- 
tics are run only once the initial linear programming relaxation is solved: 
Polynomial-time reduction strategies (Lin, 1965) use the relaxed solution 
to identify a subset of variables whose values should be fixed. Typically, 
and inevitably, the dive will remove the optimal solutions; the optimum of 
the reduced instance is not necessarily optimal in the original instance. It 
is hoped, nevertheless, that improving integer feasible solutions can still be 
extracted from solutions of the modified instance. Despite substantial re- 
cent progress, we believe that there is still considerable scope for further 
work in this general direction. 

The approach we propose is based on the observation that there is a com- 
plex structure to many real-world problems. The objective function is usu- 
ally highly structured in itself, and hence not best thought of as the generic 
c T x form, which is the usual starting point of theoretical analyses. In many 
instances, the overall objective function is a linear combination of many 
disparate terms with each term representing a total penalty for some un- 
wanted structure within a possibly small part of a solution, or profit from 
some desirable structure. Furthermore, the different terms in the objective 
are rarely equally hard to handle. Some can be relatively easy, but some 
terms can be troublesome; representing them can require a large number of 
auxiliary variables. For example, some term might well contain a penalty 
that uses a "sum of absolute values" or "max of a max" structure that can 
require many extra variables to compute. Hence, the first ingredient of our 
approach is to generate sub-instances by selectively keeping some terms 
in the objective. Which terms should be kept in the objective, will become 
clearer once we sketch out the rest of our approach. 
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The second ingredient is the aggregation of variables: Defining a new vari- 
able to be some linear or simple non-linear function of existing variables, 
such as the maximum of a small subset. In many cases, a human expert 
inspecting the sub-instance would suggest that particular aggregations of 
variables would allow it to be rewritten into a form that is much easier 
to solve. Current diving heuristics, however, do not perform any aggrega- 
tions that are specific to the diving, despite the potential to greatly reduce 
the number of variables and, furthermore, to do so in a manner that is quite 
different than reductions by means of fixing values of variables. We distin- 
guish different potential intended usages of aggregation as follows: 

• Equivalence-based. The new variable is 'logically equivalent to the old 
variables' in the sense that it can be used to replace them, and still pre- 
serve optimality and validity. This is the type of aggregation that is per- 
formed in preprocessing. 

• Validity-based. The new aggregate variable is not intended to replace the 
old variables fully: They either supplement the original ones, or can be 
used to provide valid cuts and valid lower bounds. 

• Solution-based. The new aggregate variable is not intended to replace 
the old variables fully: They either supplement the original ones, or can 
be used to provide valid solutions (upper bounds). This, in a certain 
sense, can be thought of as the dual of validity-based aggregation. 

This concept of validity-based aggregation underlies many of the construc- 
tions in this paper. 

Finally, current diving heuristics seem to have linear programming relax- 
ation as the input and (possibly) an integer feasible solution as the output. 
We suggest a different approach, (possibly) producing integer feasible solu- 
tions together with lower bounds for the original instance, without utilising 
its linear programming relaxation. This makes it possible to run the div- 
ing heuristic prior to computing the first linear programming relaxation of 
the original instance and improve the performance of subsequent cut gen- 
eration. In step one of our approach, a much smaller ("surface") instance 
of integer programming is derived by picking suitable components of the 
original instance and applying the validity-based aggregation introduced 
above. By solving this instance, we obtain a valid lower bound for the orig- 
inal instance, in addition to a feasible input to the second step. In step two, 
we dive into a neighbourhood obtained by reducing the feasible solution of 
step one. Notice that the ("dive") neighbourhood can often be much smaller 
compared to the original search space, whilst still being guaranteed to con- 
tain a feasible solution to the original instance. Various strategies can be 
used to control the order of performing the individual dives. In the long 
run, it would make sense to explore the possibilities of applying various 
strategies to initialise and perform dives, automating the choice of an aggre- 



3 



Offers / Model 


Monolithic 


Surface 


Dive 


Lower bounds 


/ 


/ 




Upper bounds 


/ 




/ 


Completeness 


/ 







Fig. 1. A schematic overview of the suggested decomposition: As a primal heuristic 
for the convergent monolithic model, we employ an aggregated surface model and 
a series of dives, bounded by the solution obtained at surface. 

gating reduction strategy, as well as of tailoring automated reformulations 
better to instances obtained in such reduction strategies. 

1.2 The Result 

In this paper, we make an initial step in this direction, giving a "proof of 
concept" implementation for a particular class of instances from timetabling. 
In step one of our approach, a much smaller instance of integer program- 
ming is solved. Feasible solutions of the smaller instance define the neigh- 
bourhoods to dive into in step two. The smaller instance is derived us- 
ing the "validity-based aggregation", introduced above. We also hint at the 
power of so-far-non-automated reformulations on instances obtained using 
such reduction strategies. 

In particular, we propose such a heuristic decomposition for Udine Course 
Timetabling, a benchmark course timetabling problem, used in Track 3 of 
the International Timetabling Competition 2007 (Di Gaspero and Schaerf, 
2006). There, as in many other timetabling problems, feasible solutions are 
determined by a bounded graph colouring component (Welsh and Powell, 
1967; Burke et al., 2004), more specifically by an extension of a pre-colouring 
bounded in the number of uses of a colour, which is difficult in its own right 
(Even et al., 1976; Garey and Johnson, 1979; Bodlaender and Fomin, 2005). 
The quality of feasible solutions is measured by the number of violations of 
four additional complex soft constraints. These soft constraints place em- 
phasis on: 

• suitability of rooms with respect to their capacities 

• suitability of the spread of the events of a course within the weekly timetable 

• minimisation of the number of distinct rooms each course uses 

• desirability of various patterns in distinct individual daily timetables. 

Note that the latter three soft constraints will be formulated using the max 
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function. Modelling of the soft constraints, as opposed to the bounded graph 
colouring component on its own, entails an increase in the dimensions of 
the model, and consequently of the run-time, making even modest real- 
life instances difficult to solve using a stand-alone integer programming 
solver (Burke et al., 2008b, 2007), or branch-and-cut procedures (Avella and 
Vasil'ev, 2005; Burke et al, 2008a). 

Hence, we propose a heuristic that decomposes the process of solving the 
instance into two stages, which are outlined in Table 1. The "surface" stage 
is restricted to the bounded colouring problem with the two period-related 
soft constraints. The "dives" fix certain features of the solution found at the 
surface, and produce a locally optimal solution to the complete instance. We 
employ two "control strategies" for executing the dives. In "anytime strate- 
gies", dives are executed whenever a feasible bounded colouring is found 
at the surface. In "contract strategies", a time limit for the search at the sur- 
face is given, and dives are executed only afterwards. As it seems, previ- 
ously underutilised "contract strategies" with multiple "neighbourhoods" 
searched in the dives in a suitable order perform considerably better than 
the traditional anytime strategies using a single type of restriction. 

A solver based on this approach was implemented using ILOG Concert and 
ILOG CPLEX 11. Using non-automatically reduced integer programming 
formulations, it yields good solutions to instances of up to 434 events and 
81 distinct enrolments ("curricula") within thirty minutes of run time on 
a desktop PC. The lower bounds obtained at the surface are better than 
those which Lach and Liibbecke (2008a) obtain within the same time limits. 
Overall, the heuristic produces solutions for the Udine Course Timetabling 
Problem, together with bounds on their distance from optimality. 

The two-fold aims of this paper are reflected in its structure. The immedi- 
ate goal is to obtain better solutions and lower bounds to the Udine Course 
Timetabling problem. The problem is introduced in Section 2. The particu- 
lar heuristic is outlined in Section 4. The integer programming formulations 
of the subproblems we employ are presented in Section 5. Finally, computa- 
tional experience is described in Section 7. We regard this as part of a longer- 
term goal, however, which is the development of methods to better control 
and exploit various, possibly automated, decompositions and reformula- 
tions. The methodology of "Multiphase Exploitation of Multiple Objective- 
/value-restricted Submodels" (or "MEMOS" for short) is outlined in Sec- 
tion 3. We are not claiming that any of the individual techniques presented 
in Section 3 are novel per se. Rather, we are providing a basis for a rationali- 
sation and classification of previously known methods, which are surveyed 
in Section 6. Both the particular solver and the more generally applicable 
approach to the design of hybrid metaheuristics for complex problems may 
be of separate interest. 
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2 Problem Description 



2.1 University Course Timetabling 

In general, timetabling problems share the search for feasible colourings of 
conflict graphs, where vertices represent events and there is an edge be- 
tween two vertices if the corresponding events cannot take place at the 
same time (Burke et al., 2004). This graph colouring component is often 
bounded in the number of uses of a colour and has to provide an exten- 
sion of pre-coloured conflict graph, such that the total number of violations 
of certain soft constraints is minimised. In university course timetabling 
(Bardadym, 1996; Burke et al, 1997; Petrovic and Burke, 2004; McCollum, 
2007), these soft constraints often stipulate that events should be timetabled 
for rooms of appropriate sizes. At least or at most a certain number of days 
of instructions should be timetabled for groups of students and individual 
teachers, and daily timetables of students or teachers should not exhibit 
particular patterns. For example, a single event per day or long gaps in 
a daily timetable, or on the other hand, six events per day with no gap 
around lunch time may be deemed undesirable. The particulars of each 
problem instance vary widely from university to university, and a num- 
ber of both exact and heuristic search methods have been attempted on 
a range of instances. For surveys, see Carter and Laporte (1998); Schaerf 
(1999); Burke and Petrovic (2002). Out of the numerous variations of the 
problem in use, instances from the University of Udine (Di Gaspero and 
Schaerf, 2003) and Purdue University (Rudova and Murray, 2003; Murray 
et al., 2007), together with random problem instances used in the Interna- 
tional Timetabling Competition 1 in 2002 and 2007 have recently started to 
be considered as benchmark problems. 

2.2 The Udine Course Timetabling Problem 

The Udine Course Timetabling problem studied in this paper is maintained 
by Di Gaspero and Schaerf (2003) at Universita degli studi di Udine. There are 
two important assumptions: 

• Events are partitioned into disjoint subsets, called courses; events of any 
one course have to take place at different times, are attended by the same 
number of students, and are freely interchangeable 

1 See Di Gaspero et al. (2007) or http://www.cs.qub.ac.uk/itc2007/ 
for International Timetabling Competition (ITC) 2007 and 

http: / /www. idsia. ch/Files/ttcomp2 02/ for ITC 2002. 
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Fig. 2. The colouring component of a trivial instance of course timetabling: There 
are three courses, Juggling, MathlOl and AlgolOl, with 1, 4 and 3 events respec- 
tively. One group of students attends both MathlOl and AlgolOl, another group 
attends both Juggling and AlgolOl, and events attended by a group of students 
cannot overlap in time. Each vertex in the associated conflict graph represents an 
event and there is an edge between two vertices, if the two corresponding events 
cannot take place concurrently. The edge is dashed, if this requirement is given by 
the enrolment of a student in both courses. 

• A small number of distinct, possibly overlapping, sets of courses, repre- 
senting enrolments prescribed to various groups of students, are identi- 
fied and referred to as curricula. 

Due to the second assumption, this problem is often referred to as "curricu- 
lum based timetabling", as opposed to "student enrolment based timetabling", 
which tries to minimise the number of conflicts among a possibly large 
number of enrolments. See Figure 2 for an illustrative example. The com- 
plete input can be captured by seven constant sets and eight mappings: 

• C, U, T, R, D, P are sets of courses, curricula, teachers, rooms, days, and 
periods, respectively 

• U is the non-empty set of courses in curriculum u 

• N is a subset of C x P, giving forbidden course-period combinations 

• E c is the number of events course c has in a week 

• S° is the number of students enrolled in course c 

• Q c is the prescribed minimum number of distinct week-days of instruc- 
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tion for course c 

• T is the subset of courses C taught by teacher t 

• A r is the capacity of room r 

• L> d is the subtuple (ordered subset) of P corresponding to periods in day 

d 

• W is a vector of non-negative weights (W RCap , W TSpr , W TCom , W RStb ) for 
the four soft constraints described below. 

Informally, the goal is to produce a mapping from events to period-rooms 
pairs such that: 

(1) For each course c, E c events are timetabled 

(2) No two events take place in the same room in the same period 

(3) No two events of a single course, no two events taught by a single 
teacher, and no two events included in a single curriculum are taught 

Table 1 

Instances of the Udine Course Timetabling problem: numbers of rooms and peri- 
ods; the number of courses and the sum of their events in a week; frequency, or the 
portion of period-room slots in use, and utilisation, or the portion of period-seat 
slots in use (Beyrouthy et al v 2007); the number of distinct enrolments ("curric- 
ula"); numbers of edges and density in conflict graphs (CG) with vertices repre- 
senting courses, rather than events (Burke et al., 2007). 
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The objective is to minimise a weighted sum of penalty terms 
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where: 



• p RCa P (for "room capacity") is the number of students left without 
a seat at an event, summed across all events; this is the value of the 
number of students attending an event minus the capacity of the al- 
located room over all events where the value is positive 

• p TS P r (for "spread of events of a course over distinct week-days") 
sums the value of the number of prescribed distinct week-days of in- 
struction minus the actual number of distinct week-days of instruc- 
tion over all courses where the value is positive 

• p TCom (for "time compactness") is the number of isolated events in 
daily timetables of individual curricula; "[f]or a given curriculum we 
account for a violation every time there is one lecture not adjacent to 
any other lecture on the same day" (Di Gaspero et al., 2007) 

• P RStb (for "room stability") is the number of distinct course-room al- 
locations on the top of a single course-room allocation per course. 

In the original paper of Di Gaspero and Schaerf (2003), there were described 
only four instances of the Udine Course Timetabling Problem of up to 252 
events and 57 distinct enrolments, with (W RCap , W TSpr , W TCom , W RStb ) = 
(1, 5, 2, 0). Fourteen more instances of up to 434 events and 81 distinct en- 
rolments have now been made available in Track 3 of the International 
Timetabling Competition, with weights (W RCap , W TSpr , W TCom , W RStb ) = 
(1, 5, 2, 1). Their dimensions are summarised in Table 1. 

2.3 An Integer Programming Formulation 

Formally, the Udine Course Timetabling problem can be described using 
an integer programming model. Doing so necessitates the choice of deci- 
sion variables. This is of little importance if the sole purpose is to formu- 
late the problem formally, but becomes of paramount importance, if the 
performance of the formulation is evaluated. It seems tempting to see the 
problem as a variation of the three-index assignment and to use binary vari- 
ables for each event-room-period combination. This corresponds to the triv- 
ial formulation of graph colouring, where the number of binary variables 
is the product of the number of vertices and an upper bound on the num- 
ber of colours. There are, however, many alternative formulations of graph 
colouring (Burke et al., 2007) and it seems reasonable to use the best avail- 
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able formulation of graph colouring that would admit formulation of the 
soft constraints. After exploring a number of such alternatives, Burke et al. 
(2007) proposed a formulation based on a suitable clique-partition, which 
is given implicitly in many graph colouring applications. For Udine Course 
Timetabling, this formulation translates to a smaller number of "core" bi- 
nary decision variables x, given by the product of the numbers periods, 
rooms, and courses, rather than events. Outside of those, there are depen- 
dent variables v, q, z, and y, whose values are derived from the values of x 
in the process of solving: 



• x are binary decision variables indexed with periods, rooms and courses. 
Their values are constrained so that in any feasible solution, course c 
should be taught in room r at period p, if and only if x p>rtC is set to one. 

• v are binary decision variables indexed with courses and days. Their val- 
ues are constrained so that in any feasible solution, there is at least one 
event of course c held on day d, if and only if v dyC is set to one. 

• q are integer decision variables indexed with courses, whose values are 
bounded below by zero and above by the number of days in a week. Their 
values are constrained so that in any feasible solution, q c is the number of 
days course c is short of the recommended days of instruction, Q c . 

• z are binary decision variables indexed with curricula, days, and an index- 
set Z of natural pattern-penalising constraints. Their values are constrained 
so that in any feasible solution, z UjdjS is set to one if and only if the pattern- 
penalising constraint indexed by s e Z is violated in the timetable for 
curriculum u and day d given by the solution. 

• y are binary decision variables indexed with rooms and courses. Their 
values are constrained so that in any feasible solution, y TjC is set to one if 
and only if room r is used by course c. 



The objective function can be expressed as: 



min W RCa P EE E ^,,c (S c - A') + w TCom E E E *M, 
reRpeP cec ueUdeDsez 

S c >A r 

+ w TS - E*=+ w RStb E { { E yr, c ) - A 

cec cec V VreR / / 
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Hard constraints can be formulated as follows: 
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Constraint (1) enforces a given number of events to be taught for each 
course. Constraint (2) ensures no two events are taught in a single room 
at a single period. Constraints (3-5) stipulate that only one event within a 
single course or curriculum, or taught by a single teacher can be held at any 
given period. Notice the similarity of constraints (1-5) and the clique-based 
formulation of bounded graph colouring (Burke et al v 2007). Finally, con- 
straint (6) forbids the use of some periods in timetables of some courses, 
corresponding to a pre-colouring extension. Notice also that constraint (5) 
renders constraint (3) redundant, if there are no courses outside of any cur- 
ricula. 



The formulation of soft constraints is less trivial and, as will be shown in 
Section 7 and Table 4, the formulation presented below proves to be quite 
challenging for modern general purpose integer programming solvers, even 
after the strengthening proposed by Burke et al. (2008a). Values of v are 
forced to maxima of certain subsets of x using constraints (7-8), in effect 
constructing daily timetables for individual curricula. Subsequently, the 
number of distinct week-days of instruction that course c is short of the 
prescribed value Q c , can be forced into q c using constraint (9): 

VcG CWGDVpGL> d J2 x P,r,c <v d ,c (7) 

reR 

VcGCVdGD E X P,r,c >Vd,c (8) 

reR P eD d 

VcgC J2 v d,c >Q c -q c (9) 

deD 



The "natural" formulation of penalisation of patterns (Burke et al., 2008b) 
occurring in daily timetables of individual curricula goes through the daily 
timetables bit by bit, first by checking isolated events in the first and the 
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last period of the day, and later looking for triples of consecutive periods 
with only the middle period occupied by an event: For an instance with 
four periods per day this is: 

\/u G U,tZ G D,V(pi,p 2 ,P3,P4) £ D d J2 Y ( X pi,r,c - Xp 2 ,r,c) < Z uA1 

ceu u reR 

(10) 

Wu G U,d G B,W (p 1 ,p 2 ,P3,P4) e D (Xpt,r,c - Xpa,r,c) < Zu,d,2 

ceu lt reR 

(11) 

Vlt G U,d G D,V(pi,j9 2 ,P3,P4) e -D ( X P2,r, C - X Pl ,r,c ~ Xp 3 ,r,c) < Z U ,d,3 

ceu u reR 

(12) 

Vlt G U,d G D,V(pi,j9 2 ,P3,P4) G -D J! 51 0W,c - ^p 2 ,r, C ~ ^W.c) < Z Uj(J ,4 

(13) 

Just as for constraints (7-9), cuts strengthening constraints (10-13) were de- 
scribed by Burke et al. (2008a). Finally, values of y, similar to values of v, 
are forced to maxima of certain subsets of x: 

VpGPVrGRVcGC x p ^ c <y r>c (14) 

Vr G RVc G C x p,r,c > Vr,c (15) 

These constraints complete the formulation, which will be later referred to 
as Monolithic. 



2.4 How Difficult is this for Exact Solvers? 



One might hope that it is enough to pass this formulation to a modern in- 
teger programming solver and wait. After all, integer programming has 
been used in timetabling ever since Lawrie (1969) generated feasible solu- 
tions to a school timetabling problem using a branch and bound procedure 
with Gomory cuts. Tripathy (1984) and Carter (1989) solved instances of 
a course timetabling problem of up to 287 events with several soft con- 
straints using Lagrangian relaxation. They have not, however, introduced 
any constraints penalising interaction between events in timetables other 
than, of course, straightforward conflicts. More recently, a number of mod- 
est instances of course timetabling problems have been tackled using off- 
the-shelf solvers, without introducing any new cuts (Dimopoulou and Mil- 
iotis, 2004; Qualizza and Serafini, 2005; Daskalaki et al., 2004, 2005; Mirhas- 
sani, 2006). For instance, Daskalaki et al. (2004, 2005) solved instances of 
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up to 211 events using ILOG CPLEX. Of special interest are studies by 
Al-Yakoob and Sherali (2007) and Schimmelpfeng and Helber (2007), who 
have modelled a larger number of constraints. In the most rigorous study 
so far, Avella and Vasil'ev (2005) presented a branch-and-cut solver for the 
Benevento Course Timetabling Problem, which forbids some interactions 
of events in timetables other than conflicts using hard constraints. In this 
setting, Avella and Vasil'ev (2005) have been able to solve instances of up to 
233 events and 14 distinct enrolments, but conceded that application of their 
solver to the four small instances of Udine Course Timetabling available in 
2005 yielded "poor results". In Udine Course Timetabling, such interactions 
are penalised by soft constraints (7-15), which (it turns out) make the prob- 
lem considerably more difficult. Several integer programming formulations 
of the problem have been studied by Burke et al. (2007, 2008b), and more 
recently by Lach and Liibbecke (2008a,b). Both approaches can now solve 
the original instances of Udine Course Timetabling from 2005, but optimal 
solutions to instances from the International Timetabling Competition 2007, 
or optima for large real-life instances of university course timetabling prob- 
lems are, as it seems, out of reach, so far. 

In particular, it seems that it is the run-time of the linear programming (LP) 
solver that is prohibitive to progress on more difficult instances from the 
International Timetabling Competition (ITC) using the Monolithic formu- 
lation. For relatively easier instances (compOl and compll from ITC 2007), 
root relaxation obtained from Monolithic takes less than thirty seconds to 
solve using the default ILOG CPLEX 11 Dual Simplex LP Solver and the 
search proceeds swiftly. On the remaining instances, however, the root re- 
laxation is far more expensive, taking up to 6489 seconds on comp07. For 
complete results obtained using the CPLEX 11 Dual Simplex LP Solver, see 
Table 4. Informal evidence suggests that run times of both the CPLEX 11 
Barrier and the CPLEX 11 Primal Simplex LP solvers are shorter, whereas 
run-times of either SoPlex (Wunderling, 1996) and CLP (Forrest et al., 2004) 
are longer. Hence, it is not particularly surprising that for most instances 
(all except from compOl, comp05, and compll), CPLEX Mixed Integer Pro- 
gramming (MIP) solver using Monolithic and default settings does not pro- 
duce any feasible solution within 40 CPU units, branching at the speed of 
ten to thirty nodes per hour. This can hardly be described as satisfactory 
progress. 

Although the total run-time or time spent at the root node may not be of 
particular importance, as long as the instance is solved close to optimal- 
ity within a reasonable time limit, for example over a weekend, relying on 
robustness of modern solvers does not seem to be an option for real-life 
instances (Murray et al., 2007). First, such instances are often many times 
larger than those used in the International Timetabling Competition 2007, 
including comp07. Second, their formulations are considerably more com- 
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plex and tend to change over time, whereas the behaviour of general solvers 
on monolithic formulations tends to be rather difficult to predict, with run- 
time often rising rapidly on addition of seemingly trivial constraints. Fi- 
nally, general solvers only seldom perform in the "anytime" fashion of Zil- 
berstein (1996) on monolithic formulations, which would make it possible 
to obtain a solution whenever the solver is stopped. It thus makes sense 
to study decompositions of the problem and the performance of heuristic 
methods based on them. 



3 Multiphase Exploitation of Multiple Objective-AValue-restricted Sub- 
models (MEMOS) 

In short, our approach can be seen as a hybridisation of loosely coupled in- 
teger programming solvers working on very different sub-problems of sim- 
ilar difficulty. Each solver is, time-limit permitting, exact within the given 
search space, and together, the solvers provide heuristic solutions with global 
lower bounds. It seems, however, useful to set this description in the broader 
context of hybrid metaheuristics, even if this still seems difficult to do pre- 
cisely, but concisely, despite the numerous attempts (Calegari et al., 1999; 
Puchinger and Raidl, 2005; El-Abd and Kamel, 2005; Raidl, 2006) to develop 
a taxonomy of the field. 

3.1 Hybridisation 

Hybridisation can clearly occur in two different fashions: 

• Algorithm hybridisation, using multiple algorithms for the same prob- 
lem 

• Formulation hybridisation, using multiple subproblems or multiple lev- 
els of abstraction 

Although, in theory, hybridisation can occur in both fashions, in practice, 
heuristics often exploit only algorithm hybridisation: Multiple algorithms 
producing solutions to the complete problem, with or without a feasible 
solution on the input. A typical example of algorithm hybridisation is the 
"construct-improve" scheme, which is a standard practice in metaheuristic 
solvers. A constructive algorithm is used in order to produce one, or many, 
different initial solutions, which are then improved using some variant of 
local search. See the concise overview by Blum and Roli (2003) for details. 
Embedding of multiple primal and improvement heuristics (Fischetti and 
Lodi, 2003; Danna et al., 2005; Berthold, 2007b) in general integer program- 
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Fig. 3. A schematic illustration of the data flows in the metaheuristic. 



ming solvers can be seen as another example of algorithm hybridisation. In 
contrast, the deliberate creation and exploitation of multiple sub-problems, 
restricted to certain components in the objective function, with certain val- 
ues fixed, or both, seems underutilised. 



3.2 Submodels 



The critical point in our approach is precisely that: sub-problems solved 
at various stages are different, although possibly overlapping. That is, we 
create sub-problems using combinations of 

• Value-Restrictions: Fix the values of a subset of the variables, distance 
from a solution with respect to Hamming distance (L norm), or sums 
across subsets of variables. This is a standard method to generate a neigh- 
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bourhood. 

• Objective-Restrictions: Suppose that we have a minimisation problem 
and some penalty term P R is part of the objective, with associated weight 
W R . We then use one of two options: 

• Objective-Ignore-Term We drop all consideration of the penalty - this 
corresponds to setting W R = 0. 

• Objective-Fix-Term We force P R = 0. This converts the penalty to a 
hard constraint. 

The objective restriction might loosely be thought of as the dual of the 
variable restrictions. 



In any case, the intent is that the restriction will give a sub-problem that can 
be significantly simplified compared to the full problem, and of comparable 
difficulty to other problems obtained in the decomposition. 

This can be illustrated on a class of problems, which often occurs in Schedul- 
ing and Timetabling applications, where there is a graph colouring compo- 
nent expanded over another dimension. Imagine for example the assign- 
ment of jobs to time-periods, expanded over machines, or the assignment 
of events to time-periods, expanded over rooms. In such cases, only one 
part of the objective is usually related to the expanded problem, whereas 
the rest is related to the graph colouring component. This results in struc- 
tures such as: such as: 



nop \ 

min CjXj + dj max z fc + y^ di max(/(Z, maxx m ), 0)1 (16) 

K i=l j=l 3 l=o-\-l ) 

maxx k < 1 Vj G {1,2,..., o} (17) 

Ax>b (18) 

xe{0,l} n (19) 



where n, m, o, and p are integral constants, x is a vector of n integral de- 
cision variables, A is an n x m constant matrix, b, c, and d are compatible 
constant vectors, and for each j = {1,2,..., p}, Ij is a subset of {1, 2, . . . , n}. 
Notice that the second two terms in the objective function, weighted by d, 
are related only to the graph colouring component, and necessitate imple- 
mentation of the max function using additional variables. The number of 
"core" variables x can be much less than the number of additional vari- 
ables implementing the max function over various subsets of x. Hence, it 
makes sense to use the restricted objective: 



I dj max Xk + y^ di max(/(Z, maxx m ), 0) , (20) 
u=i fce/j i=o+i me/i 
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ignoring the first term in the original objective (16). Given the constraint 
matrix A exhibit certain structure (17), this should allow for aggregation of 
n variables x into o variables y: 

min [jz^y, + djmaxC/tf, y fl(J) ), 0) I A'y > b' (21) 

\j=l l=o+l J 

y G {0, 1}° (22) 

where g is a suitable mapping. Notice that some constraints may have to be 
relaxed in deriving the new constraints (A'b') from the original ones (Ab). 
In Table 1 and below, we denote such a subproblem as "surface". Such an 
objective-restricted subproblem, however, should solve considerably faster 
than the original problem (16-19). 

Once a solution to the objective-restricted subproblem is found, it makes 
sense to resolve the original problem (16-19) with the additional constraints: 

maxx k = y j Vj G {1, 2, . . . , o} (23) 
m&x(f(l,maxx m ) = yi V7 G {o + 1, o + 2, . . . ,p}. (24) 

mdJi 

In Table 1 and below, we denote such subproblems as "dives". This should 
result in considerable reductions of the model in presolving, and conse- 
quently in value-restricted subproblems which solve considerably faster 
than the original problem (16-19). 

One can also use any combination of these, and could also separately treat 
different terms of the objective. This is, however, limited by the availabil- 
ity of a decomposition into such sub-problems, given by the structure of 
a solution and the objective, applicable reformulations, and our ability to 
automate their execution. For many "academic" problems there tends to 
be just a single term on the objective - for example, for the TSP we have 
only the tour length, and the objective restriction is hence not relevant. For 
many real-life problems, however, it is quite likely that the solution has a 
more complex structure, and there are multiple soft constraints correspond- 
ing to multiple penalty terms in the objective. University course timetabling 
(Murray et al., 2007) as opposed to pure graph colouring, or vehicle routing 
problems (Chabrier, 2006; Pisinger and Ropke, 2007) as opposed to the trav- 
elling salesman problem, provide good examples of possible application 
areas. In this case ignoring or fixing some terms can give significantly eas- 
ier problems, yet still provide enough information to guide solution of the 
full problem. As we will see, automated preprocessing in-built in modern 
integer programming solvers is often too weak, and we do need to apply 
non-automatic reformulations "by-hand". Ultimately, our work is directed 
towards automating decompositions and reformations of such problems 
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search at dive 1 search at the surface level dive 2 search at the surface level dive 3 run-time 
the surface (continued) (continued) 



Fig. 4. A schematic illustration of information flows in the anytime algorithm: fea- 
sible solutions at the surface are transformed into neighbourhoods for dives, exe- 
cuted immediately. Upper bounds (costs of best-so-far solutions) are being carried 
over from dive to dive. 



with "multiple soft constraints", exploiting cases in which the objectives 
can have differing levels of importance. 



3.3 Control Strategies 



Another important and often neglected question in the design of hybrid 
metaheuristics asks for the sequence of periods of searching at the surface 
and periods of diving in neighbourhoods, sometimes referred to as "the 
control strategy" (Puchinger and Raidl, 2005; Raidl, 2006). Very often, a 
variant of the anytime algorithm (Zilberstein, 1996) is used, perhaps be- 
cause it fits well into the large neighbourhood search scheme. There are, 
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however, two obvious alternatives of the producer-consumer relationship: 

• Anytime algorithm: search at the surface and dives are interleaved and 
dives are perfomed as soon as suitable feasible solutions are found. The 
algorithm is parametrised with the types of dives performed, their order, 
and possibly the onset of execution of each type. See also Figure 4. 

• Contract algorithm: search at the surface is carried out for a certain time 
and only then are the dives performed, firstly ordered by dive type with 
the least expensive dives coming first, and secondly ordered in the as- 
cending order of the cost (in minimisation problems) of the solution at 
the surface the dive is derived from. Hence, the algorithm is parametrised 
with the types of dives performed, the time limit and the proportion of 
the time limit spent at the surface. See also Figure 5. 

It seems that, in practice, a time-limit is often given and considerable per- 
formance improvements can be gained from exploiting this in a Contract 
algorithm. If the neighbourhoods can be searched in the ascending order 
of the cost of the solution at the surface and the cost of the present-best so- 
lution found in dives is applied as an upper bound in all dives except the 
first one, a considerable number of dives can be cut off immediately and 
the algorithm needs not to be parametrised with the onset of diving. (This 
is suggested in Figure 5.) Notice that the ascending order of the cost cor- 
responds to the inverse order of the time of discovery, if we are producing 
only progressively improving solutions at the surface. 

Generally, such Multiphase Exploitation of Multiple Objective- /value-restricted 
Submodels (MEMOS) will no longer be complete or exact. The intention is, 
however, to produce good solutions quickly. 



4 MEMOS for Udine Course Timetabling 

In the design of the particular MEMOS heuristic, we need to make specific 
decisions as to which sub-problems the "surface" and "dive" components 
need to solve. In doing so, it seems useful to think of the problem from the 
multi-objective perspective. In the Udine Course Timetabling problem in- 
troduced in Section 2, a given set of weights is used to linearise the penalties 
given by the number of violations of four soft constraint ("objectives") into 
a single-objective problem. The four objectives are not necessarily linked, 
however. The formulation introduced in Section 2.3 models the four soft 
constraints one-by-one, each objective using a separate subset of constraints 
and a separate subset of variables derived from the main decision variables. 
Hence, it is natural to consider sub-problems arising from setting weights 
related to a subset of the four soft constraints to zero. 
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search at search at the surface level search at the surface level dives run-time 

the surface (continued) (continued) 



Fig. 5. A schematic illustration of the contract algorithm: dives of all types are post- 
poned until the surface search finishes (runs out of time), when they are carried out 
in the descending order of the time of discovery of the neighbourhood, with the 
upper bound found from in the best neighbourhood often cutting off the dives into 
other neighbourhoods as soon as their local lower bounds are established. 



If some weights are set to zero, the corresponding objective can be ignored 
and the problem formulation can be simplified. Ideally, one would be able 
to rely upon the preprocessing in-built in an integer programming solver 
to do this. As suggested by Table 2, however, present-day general purpose 
integer programming solvers are not particularly effective for such major 
reformulations. Non-automatic reformulation can, however, result in much 
faster relaxations, outside of having other desirable properties (Trick, 2005). 
In our case, such a desirable property might be the ease of conversion of 
solutions obtained at the "surface" to neighbourhoods for the dives. As 
suggested in Table 4, a careful choice of the decomposition and manual 
reformulation can lead to the run-time of the linear programming solver at 
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Fig. 6. A schematic illustration of the penalties considered in the formula- 
tions we use: Monolithic uses the complete objective function, Surface sets 
w RCa P = w RStb = Q/ DayF i xe d has P TSpr fixed, PeriodFixed has P TSpr and p TCom 
fixed, and DayDecomp has P TS P r fixed and W RStb = 0. 



the relaxation either in the surface component or within a dive being cut by 
factor of more than two thousand, compared to the "monolithic" formula- 
tion. 



4.1 Objective-restricted Submodels 



One could consider generating reduced problems by setting any subset of 
the weights to zero. With four penalties, there are 2 4 — 1 = 15 possible non- 
trivial choices. However, we only use the following problem: 

• Surface: We set W RCap = W RStb = so that room over-filling and room 
stability are not considered. Since no explicit room assignments are re- 
quired, there can be much fewer variables. We will see later that the 
ternary variables x p r c are no longer required, which gives us much shorter 
run-times of the linear programming solver. (See Table 4.) However, we 
add extra constraints that bound the number of rooms used in any single 
period; they impose that there will always be enough rooms for courses 
at each period, and so guarantee that a surface solution can always be 
extended by the dive to a feasible solution to the full problem. However, 
these extra constraints convert the colouring to bounded colouring and 
make solving the "surface" search harder, both in theory and practice. 

The sets of the objectives considered by each formulation are illustrated in 
Figure 6. This also serves to emphasise the crucial point that the methods 
we use are intended for the case of multiple terms in the objective, or equiv- 
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alently, for multiple classes of relatively-complex soft constraints. 

One can reasonably ask why we use this particular split between surface 
and dive models. Firstly, the graph colouring component is important in 
the instances under consideration. In Table 1, we can see that the ratio of 
events to available room slots is typically 60-90% which is quite high, so it 
is reasonable that the bounded colouring problem is hard so solve, and so 
ought to be solved first. That is, the bounded colouring needs to be done 
at the surface level. Also, Table 2 gives optimal objective vectors for dif- 
ferent weight vectors. (It does this for the small instance compOl that hap- 
pens to be fairly easy to solve exactly, and uses the monolithic IP formula- 
tion detailed in the previous section.) We see that the time-related penalties 
(P TSpr and p TCom ) make the largest contributions to the final objective, and 
so it is natural to address these first. Limited experimentation has been con- 
ducted with solvers where W RCap = W RStb = W TCom = 0. As is suggested 
by Table 2, the solver at the surface is not considerably faster, when only 
automatic reductions are considered. Rather surprisingly, the solver is not 
considerably faster, however, even when non-automatic reformulations are 
applied, and these gains do not offset the need to spend more time div- 
ing. (See Section 5.2.) There are also alternatives exploiting interchangeable 
rooms. 



4.2 Exploiting Interchangeable Rooms: "Multi-rooms" 

Suppose that we have two rooms of exactly the same size. Then the only 
way that the rooms are not interchangeable is because of the room stability 
constraint. If we take a solution and swap the rooms at any single time 
period then no constraints are broken and no objectives affected, except 
potentially the room stability. Hence, if we are using a formulation, such 
as the surface, that does not care about room stability, then we can treat 
the rooms as interchangeable. Thus, we could replace the 2 standard rooms 
with a single room with multiplicity, which we call a "multi-room". 

A multi-room of size "(m,a)" will have multiplicity m, capacity a, and will 
function exactly as m separate indistinguishable rooms. That is: It can ac- 
commodate m separate events of sizes up to a each, without any overflow. 
This mild extension is straightforward, but has the advantage that we then 
have the option to replace the r index of the x p>rjC variables by an index into 
a set of multi-rooms. Any solution using multi-rooms can be translated into 
constraints on original rooms that define a neighbourhood in the full prob- 
lem. Within this neighbourhood, there will always be a feasible solution, 
and we can easily find the one minimising p 113 * 13 , the penalty for "room sta- 
bility". 



22 



Here we use multi-rooms only as a device to reduce the number of variables 
and remove symmetries that would enlarge the search space. However, it 
is also reasonable that multi-rooms could be directly useful. Specifically, in 
practice, the room stability might well not be quite so strict as the Udine 
timetabling instances demand. In some teaching facilities, there could well 
be a suite of teaching rooms that are all of the same size, with the same 
facilities, and which are located very close to each other. In such a case, it 
could well be that there really would be no need for a room stability penalty 
for mixed use. In this case, they could well be treated exactly as a single 
multi-room with appropriate multiplicity. 

Generally, only a few rooms do have exactly the same sizes. However, it is 
often the case that there are many rooms of roughly the same size. To ex- 
ploit this situation, note that if we join two rooms of different sizes into a 
single multi-room then it can affect the room capacity objective, but again it 
will not affect the hard constraints. Hence, we do allow the creation of relax- 
ations in which various sets of rooms are joined into associated multi-rooms 
with capacity the size of the largest original room. Specifically, if we select 
a set of multi-rooms with capacities {ci\i = 1 . . . |R|} then we can replace 
them with a multi-room of size (m, max Cj}. An intermedaite case would be 
to divide the rooms into two multi-rooms, according to whether they are 
smaller or not than some selected intermediate size. We will introduce such 
a neighbourhood generator, denoted by Surface2, in Section 5. 



4.3 Value-restricted Submodels 

For the dives, we use value-restricted submodels of the Monolithic model. 
Two pure value-restricted submodels are defined by restricting the times of 
courses based on the solution obtained at the surface level: 

• PeriodFixed dives: The periods of all the courses are fixed to be the same 
as obtained from the surface solution. 

• DayFixed dives: The days of all the courses are fixed, but not the explicit 
period within the day. That is, for each course the assigned "period" is 
relaxed to just give a day assignment. 

In both pure value-restricted subproblems, the spread of events of a course 
over distinct week-days is entirely fixed, and W TSpr hence becomes irrele- 
vant. PeriodFixed dives also fix time compactness penalty 'p TCom / and thus 
only room assignment remains to be found, minimising the room related 
penalties (P RCa P and p RStb ). In DayFixed dives, the compactness penalty 
pTCom can potentially be improved. A middle ground between these 
two approaches is offered by value- and objective-restricted submodels. 
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4.4 Value- and Objective-restricted Submodels 



Two obvious value- and objective-restricted submodels can be defined by 
restricting the times of courses based on the solution obtained at the surface 
level, and not considering P RStb or by forcing P RStb to zero: 

• DayDecomp dives: p RStb is not considered and the days of all the courses 
are fixed, but not the explicit period within the day That is, for each 
course, the assigned "period" is relaxed to just give a day assignment 
and the lowest possible penalty for time compactness and room capacity 
is sought. 

• DayFixedZeroP RStb dives: p RStb is forced to be zero and the days of all 
the courses are fixed, but not the explicit period within the day. That is, 
for each course, the assigned "period" is relaxed to just give a day assign- 
ment and an assignment with zero penalty for room stability is sought. 

These two value- and objective-restricted submodels have the advantage 
that they reduce the links between variables representing individual days, 
which translates into matrices closer to the block structure and improve- 
ments in the performance of the diving integer programming solver. Many 
other subproblems, most notably variants of PeriodFixed and DayFixed 
dives using stochastic ruining (Schrimpf et al., 2000; Dueck et al., 2002), are 
of course possible, but are not studied here. Some experience with imple- 
menting and testing the presented alternatives empirically is described in 
Section 7. 



4.5 Control Strategies 



Finally, we have designed both Anytime and Contract control strategies for 
Udine Course Timetabling, using PeriodFixed dives followed by DayFixed 
dives. The details of their functioning is outlined in Figure 7. Contract algo- 
rithms turned out to perform better, as expected, especially when multiple 
solutions can be found at Surface within the given time limit. The only dif- 
ficulty seems to be the choice of the partition of short time limits between 
Surface and dives in Contract algorithms. When a time limit is given, and 
is not too strict, however, there seem to be no reasons to use Anytime strate- 
gies. 
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class AnytimeStrategy implementing Strategy: 
constructor (config) : 
keep track of config 

initialise bestSoFar to an error value 
onFeasibleSolution ( ) : 

foreach type in config . getDiveTypes ( ) : 

if time > config . getAnytimeDivingOnset (type ) : 
Neighbourhood def = getNeighbourhood (type) 
update bestSoFar with 

dive (def, config . getDiveTimeLimit (type) ) 
onTimeout ( ) : 

return bestSoFar 

class ContractStrategy implementing Strategy: 
constructor (config) : 
keep track of config 

initialise a stack of neighbourhoods n[t] 
for each t in config . getDiveTypes ( ) 
onFeasibleSolution ( ) : 

foreach type in config . getDiveTypes ( ) : 

f ixDayNeighbs .push (getNeighbourhood (type) ) 
onTimeout ( ) : 

initialise bestSoFar to an error value 
foreach type in config . getDiveTypes ( ) : 
c = 

while time > config . getTimelimit ( ) && 

! n [type ] . empty () && c < config. getRunLimit (type) : 
update bestSoFar with 

dive (n [type ] . pop ( ) , config. getDiveTimeLimit (type ) ) 
c += 1 

class SampleAnytimeSolver ( ) : 
solve ( ) : 

config . setDiveTypes ( [FixPeriod, FixDay] ) 
foreach type in config . getDiveTypes ( ) : 

come up with a suitable time limit for type 

config . setDiveTimeLimit (type, limit) 
load an instance 

initialise a solver for Surface, instance, and config 
install AnytimeStrategy (config) as a solver callback 
try solving Surface 

class SampleConstractSolver ( ) : 
solve (timeLimit ) : 

config . set TimeLimit (timeLimit ) 

config . setDiveTypes ( [FixPeriod, FixDay] ) 

config . setDiveTimeLimit (FixPeriod, time Limit /50 ) 

config . setDiveTimeLimit (FixDay, timeLimit/5) 

foreach type in config . getDiveTypes ( ) : 

config . setRunLimit (type, 3 ) 
load an instance 

initialise a solver for Surface, instance, and config 
install ContractStrategy (config) as a solver callback 
try solving Surface for timeLimit/3 

Fig. 7. Pseudocode of the control strategies in use 



5 Non-Automatic Re-Formulations of the Subproblems 

As has been described in Section 3, at the "surface" level of the heuristic 
for Udine Course Timetabling, there is an objective-restricted neighbour- 
hood generator trying to find an assignment of events to periods, such that 
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at most a given number of events takes place at any given period, but dis- 
regarding other issues of assignment of rooms, such as room capacities. 
Feasible solutions of the assignment of events to periods are relaxed into 
value-restricted neighbourhoods given by the assignment of events to days 
or periods. It would obviously be possible to use the monolithic formula- 
tion with weights (W RCap W TSpr W TCom W RStb ) = (0, 5, 2, 0) at the "surface" 
level, and to rely on the automatic pre-solving routines in a modern inte- 
ger programming solver. As will be shown in Section 7, however, improve- 
ments in the run-time of the linear programming solver, which determine 
the run-time of the integer programming solver, of several orders of mag- 
nitude can be achieved by non-automatic reformulation. 

5.1 The Search for Good Neighbourhoods 

In the formulation that we denote as Surface, the "core" decision variables 
w are accompanied by dependent variables v, q, and z, introduced previ- 
ously: 

• w are binary decision variables indexed with periods and courses. Their 
values are constrained so that in any feasible solution, course c should be 
taught at period p, if and only if w p :C is set to one. 

The objective function can be expressed as: 

min W TCom E E E + W TSpr E<? C 
ueu deD sez cec 

Hard constraints can be formulated similarly to those in Monolithic: 
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Notice that constraint (28) is the only mention of rooms in this formulation 
of neighbourhood definition. It makes the number of rooms used in any pe- 
riod smaller than |R|, which corresponds to making the colouring bounded. 
It would be possible, for example, to constrain the number of large courses 
taught in any one period to be less than the number of large rooms, with 
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suitable definitions of "large". Similar constraints, however, do not seem to 
improve the quality of neighbourhoods significantly. 

By minor changes of constraints (7-9) in the monolithic formulation, the 
array q can be constrained using: 

Vc G CVd G DVp G i) d w p , c <v c4 (30) 

VcGCVdGD ]T w p , c >v C)d (31) 

pS-D 

VcGC Y. v c4 >Q c ~q c (32) 

Inequalities 30 and 31 in effect aggregate values from w into v. Inequality 32 
then calculates the number of days the instruction is short of the recom- 
mended value in Q. 

Patterns can be penalised again using the natural formulation, used already 
in constraints (10-13) in the monolithic model. For instances with four pe- 
riods per day, we obtain: 

VwGUdG DV(pi,p 2 ,P3,P4) e D ^] (w puc - w P2 , c ) < z uA i (33) 

ceu u 

Vu G U d G DV (pi,p 2 ,P3,P4> e D ( w P4,c ~ w ps,c) < z u,d,2 (34) 

\/u E\J d EB\/ (p ll p 2l p 3l p 4 ) E D ( w P2,c - «W - w P3iC ) < z uA3 (35) 

ceu u 

\/u G XJd G DV(pi,j72,p 3 ,p4) 6 E ( w ps,c - «W - w P4jC ) < z Uidi4 (36) 

ceu 

This formulation of Surface can obviously also be strengthened by addition 
of cuts, as described by Burke et al. (2008a). 



5.2 Better Neighbourhoods? 



The approach outlined above is not the only possible neighbourhood gen- 
erator, although the range of possible formulations is somewhat limited by 
the choice of decision variables. Considering further soft constraints would 
necessitate the addition of a considerable number of variables, resulting 
in very slow LP relaxations. (See Table 4 for an illustration of the effect 
and Section 7 for discussion.) Replacing the array w with an array indexed 
with courses and days would make it impossible to implement the graph 
colouring component. The only option left, hence, seems to be the consider- 
ation of only a single soft constraints, setting (W RCap W TSpr W TCom W RStb ) = 
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(0, 5, 0, 0). This enables the removal of z and constraints (33-36). It seems, 
however, that using weights W = (0, 5, 0, 0) makes it necessary to search 
larger neighbourhoods in dives (see below for DayFixed), and results in 
worse performance on both short and long runs. Alternatively, we can ex- 
ploit interchangeable rooms. 



As has been suggested in Section 4, Surface can be thought of as the ex- 
treme case of room aggregation, in which all the |R| rooms are joined into 
a single multi-room of multiplicity \R\ and capacity the size of the largest 
room. In all these instances, the resulting capacity is always larger than the 
largest event, and so the multi-room multiplicity just gives the bound on 
the number of events per time-slot. An intermedaite case would be to di- 
vide the rooms into two multi-rooms, according to whether they are smaller 
or not than some selected intermediate size. In a neighbourhood generator 
denoted Surf ace2, the threshold room size is just taken to be the median of 
all the room sizes. This keeps the number of variables down, but still adds 
some restrictions so that larger events are likely to be evenly spread out 
over the periods. 



More formally, given the set M of pairs (m, a), representing multiplicity 
and capacity of a multi-room, the formulation we denote Surface2 uses 
the "core" decision variables x', which differ from x only in that they are 
indexed with M instead of R, together with dependent variables v, q, and 
z, introduced previously, and together with y', which are similar to y, only 
indexed with M instead of R. The objective function can be expressed as: 



mm 



in w RC * E E E <M,c (s c - a) + w TCom E E E *M, 

(m,a)GMpeP ceC weUdeDseZ 
S c >a 

+ w ts -e? c + w RStb E((E< c )-i) 

c€C cGC \ VreM / / 



Hard constraints can also be formulated similarly to those in Monolithic 



28 



(1-5): 



Vc G C 


EE 

peP rgM 


•Ep,r,c 


= E C 


(37) 


Vp G P V(m, a) G M 


E 

cGC 


•Ep,r,c 


< m 


(38) 


Vp G PVc G C 


E 

reM 


•Ep,r,c 


< 1 


(39) 


Vp G PV£ G T 


EE 


%p,r,c 


< 1 


(40) 


Vp G PVw G U 


E E 


%p,r,c 


< 1 


(41) 


V(c,p) G N 


E 


■Ep,r,c 


= 


(42) 



The remaining constraints of Monolithic can be used in Surf ace2, when x' 
replaces x, y' replaces y, and M replaces R. 

5.3 The Diving 

When an integer feasible solution is found at Surface, it waits to be trans- 
lated to a solution of the Udine Course Timetabling Problem. There are at 
least two obvious restrictions of the monolithic formulation, which we refer 
to as PeriodFixed and DayFixed dives in Section 3. It is clearly possible to 
constraint the monolithic formulation into a PeriodFixed dive by the addi- 
tion of constraints: 

Vc G CVp G P X P,r,c = Wp,c (43) 

rgR 

On the instances from the International Timetabling Competition 2007, it is 
often possible to find optima for PeriodFixed dives within seconds. As will 
be described in Section 7, however, the quality of solutions obtained this 
way varies widely from neighbourhood to neighbourhood. 

Alternatively, it is possible to dive into neighbourhoods defined by assign- 
ment of courses to days. In order to do so, values of the array D indexed 
with courses and days are pre-computed outside of the solver as: 

Vp G PVrf G D D d ' c = J2 w p,c 

P eD d 

Then, it is possible to constraint the monolithic formulation using 

VceCVdeDVpe^ ]T x p ^ c = D d > c (44) 

rGR 
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Presolving routines of modern general purpose integer programming solvers 
are typically able to shrink integer programming instances obtained in Day- 
Fixed dives by the factor of ten or more. Subsequently, optima within such 
neighbourhoods are found within minutes for the instances from the Inter- 
national Timetabling Competition. For details, see Section 7. 

5.4 Why Not Better Diving? 

A middle ground between these two solutions is represented by "ruining" 
(Schrimpf et al., 2000) into the assignment of events to days only those el- 
ements w, which seem to contribute to the objective function, or to its fac- 
tor p TCom . The first option is easy to implement, as modern integer pro- 
gramming solvers expose the pseudo-cost of individual variables, but it 
does not seem to be particularly useful for instances from the International 
Timetabling Competition 2007. The implementation of the second option 
involves a certain computational effort, but it might prove useful in some 
instances. Both kinds of "alternative dives" are disregarded in the remain- 
der of the paper, and might be the focus of future work. 
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Table 2 

The effects of disabling various subsets of soft constraints by setting the corre- 
sponding weights in the monolithic model to zero on CPLEX 10 run-time for in- 
stance compOl. Only the default in-built pre-processing and cuts are used and only 
the best feasible solution obtained within one hour per run is presented. Notice that 
when a weight is set to zero, the associated penalty is typically high. 



X 

Q 


a 

CO 

H 


2 
c 
U 
R 

£ 


cn 


LP Matrix 


Barrier LP 
Runtime 


Node with 
last feas. sol. 


Default IP 
Runtime 


pRCap 


pTSpr 


pTCom 


pRStb 


Obj. 














6639 x 5407 


44.08 s 





66 s 


2294 


30 


350 


93 














1 


6639 x 5407 


43.66 s 


370 


1066 s 


2706 


30 


350 














2 





7059 x 5757 


61.24 s 


140 


884 s 


2107 


30 





77 








5 








6665 x 5435 


40.44 s 





148 s 


2387 





350 


92 





1 











6639 x 5407 


21.68 s 





62 s 


4 


30 


350 


44 


4 








2 


1 


7059 x 5757 


67.73 s 


617 


1033 s 


2874 


30 














5 





1 


6665 x 5435 


39.41 s 


402 


1554 s 


2023 





350 











5 


2 





7085 x 5785 


66.18 s 


90 


881 s 


2365 








78 





1 








1 


6639 x 5407 


9.97 s 


180 


521 s 


4 


30 


350 


1 


5 


1 





2 





7059 x 5757 


21.21 s 


124 


518 s 


4 


30 





31 


4 


1 


5 








6665 x 5435 


10.01 s 





52 s 


4 





350 


44 


4 





5 


2 


1 


7085 x 5785 


39.50 s 


1964 


1050 s 


1672 














1 





2 


1 


7059 x 5757 


21.40 s 


487 


1416 s 


4 


30 





1 


5 


1 


5 





1 


6665 x 5435 


7.96 s 


269 


651 s 


4 





350 


1 


5 


1 


5 


2 





7085 x 5785 


22.35 s 


145 


395 s 


4 








43 


4 



Table 3 

The effects of varying the weight W TSpr on CPLEX 10 run-time for the monolithic 
model and instance compOl. 



% 

Q 
<& 
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— 

a 

CO 

R 


£ 

o 
O 

:- 

£ 


X 
CO 

£ 


LP Matrix 


Barrier LP 
Runtime 


Node with 
last feas. sol. 


Default IP 
Runtime 


pRCap 


pTSpr 


pTCom 


pRStb 


Obj. 


1 


2 


2 


1 


7085 x 5785 


19.19 s 


1053 


1205 s 


4 








1 


1 


1 


4 


2 


1 


7085 x 5785 


19.19 s 


405 


> 3600 s 


4 








2 


(6) 


1 


8 


2 


1 


7085 x 5785 


19.49 s 


1420 


> 3600 s 


5 








1 


(6) 


1 


16 


2 


1 


7085 x 5785 


19.23 s 


4822 


> 3600 s 


4 








2 


(6) 


1 


32 


2 


1 


7085 x 5785 


19.06 s 


586 


> 3600 s 


4 








2 


(6) 


1 


64 


2 


1 


7085 x 5785 


19.39 s 


3803 


2387 s 


4 








1 


5 
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6 Related Research Issues: A Brief Discussion 



It should be noted that the components of this approach draw upon a rich 
history of work in this area. Decompositions have been used in integer pro- 
gramming for a long time (Ralphs and Galati, 2006, provide a survey) and 
the effects of non-automatic reformulations are also well known (Williams, 
1978; Barnhart et al., 1993, for example). In the timetabling community, 
the "times first, rooms second" decomposition is a standard procedure. 
See Burke and Newall (1999) for references. Decompositions using integer 
programming are less common, although Lawrie (1969), for instance, used 
a form of column generation as early as 1969. Decompositions providing 
lower bounds are yet less common, and often limited in constraints they 
allow. Carter (1989) used only variables indexed with rooms and courses, 
which prohibits formulation of period-related soft constraints, but imposes 
a variant of room stability via cuts. More recently, Daskalaki et al. (2005) 
used "times first, rooms second" with integer programming models solved 
at both stages, although, in our terms, they use only a single type of dives, 
anytime control strategy, and evaluate the performance only on two small 
instances. Yet more recently, Lach and Liibbecke (2008a) have also devel- 
oped a similar approach, independently of the present authors. 2 They use 
a rather different and interesting surface component, based on their studies 
of the Partial Transversal Polytope (Lach and Liibbecke, 2008b). Their very 
respectable results are discussed below, in Section 7, and are included in 
Figure 8 for comparison. In machine scheduling, such approaches are often 
termed "machine aggregation heuristics" Leung et al. (2004) and it is not too 
difficult to translate rooms to machines and vice versa. We have, however, 
been unable to find a more general two-staged decomposition methodol- 
ogy using integer programming at both stages, comparable to the one pre- 
sented here. Overall, we heartily recommend both papers on decomposi- 
tions of timetabling (Daskalaki et al., 2005; Lach and Liibbecke, 2008a) to 
the reader, for comparison. 

Multiphase exploitation of multiple objective- /value-restricted submodels 
in some aspects resembles very large neighbourhood search (Ahuja et al., 
2002) and adaptive large neighbourhood search (Ropke and Pisinger, 2006): 
We use a number of different subproblems, which can be very large and 
slow to solve. Multiple and/ or large neighbourhoods have recently been 
very popular in local search solvers for timetabling applications. Di Gaspero 
and Schaerf (2003) and Miiller (2008), for instance, use mutiple polynomial- 



2 Results of Lach and Liibbecke (2008a) were first submitted in February 2008 
to Practice and Theory of Automated Timetabling, a conference held in Montreal, 
Canada, in August 2008. The approach of the present authors was first presented 
at the Mixed Integer Programming Workshop in August 2007. 
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time searchable neighbourhoods for Udine Course Timetabling. Abdullah 
et al. (2006, 2007) and Meyers and Orlin (2006) have proposed large neigh- 
bourhoods for timetabling, without using integer programming to search 
the them. Daskalaki et al. (2005) searches large neighbourhoods using inte- 
ger programming. Avella et al. (2007) study high-school timetabling prob- 
lem and do a form of "value-restricted diving" based on fixing the schedule 
for all but two teachers. Depending on the particular problem and imple- 
mentation, however, the "surface" component in our approach tends to pro- 
duce only few progressively improving solutions, and good lower bounds 
valid for the full problem, which is often not the case in general large neigh- 
bourhood search. In contrast to large neighbourhood search, the "dives" 
also do not have to be explored at the time of their discovery or even in 
the order of discovery, or even on the same machine. If it is possible to 
produce a number of progressively improving solutions at the surface, it 
indeed seems natural to start diving in neighbourhoods from the best so- 
lution found at the surface, with the hope that the solutions found within 
the dive can be used to cut off exploration of other dives once their lower 
bounds are known. Very large neighbourhood search nevertheless remains 
the metaheuristic approach closest to the presented one. 

There are also some similarities with the concept of "ruin-and-recreate" 
(Schrimpf et al., 2000; Dueck et al., 2002). As Lin (1965) pointed out, the key 
decision is which elements to preserve and which to relax, when the solu- 
tion is passed from the surface to a dive. This might be termed an option 
to ruin part of the solution before diving. However, it should be empha- 
sised that the overall algorithm is very different from "ruin-and-recreate". 
In particular there is no overall loop in which solutions are recycled. A dis- 
tinguishing feature or our approach is that the flow of solutions, as given in 
Figure 3, is different from the usual cycle of solutions. More specifically, in 
Figure 3, there is no arrow returning from the dive to the surface level. The 
surface merely continues to produce its next solution, without reusing the 
result of the dive directly. If the present-best solution is found in a dive, its 
value only provides an upper bound (cut-off limit) for further dives. This 
strategy seems much more compatible with the architecture of integer pro- 
gramming solvers. 

There are also connections to column generation and dimensionality reduc- 
tion in multiobjective optimisation. The heuristic surface-dive decomposi- 
tion can be likened to column generation and other exact decompositions 
studied in integer programming (Ralphs and Galati, 2006). The "surface" 
and "dive" components could also be thought of as an biobjective inte- 
ger program. Although there has been recently some progress in the study 
of biobjective integer programs (Ralphs et al., 2006), they do not seem to 
yield direct performance improvements. It would also be most interest- 
ing to study the links to dimensionality reduction (Tenenbaum et al., 2000). 
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These connections, however, are largely yet to be established, which might 
take some time: As has been pointed out by Gandibleux and Ehrgott (2005), 
the progress in understanding the true multiobjective nature of many opti- 
misation problems and related solution methods remains painfully slow. 
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Table 4 

Linear programming relaxations for instances of the Udine Course Timetabling problem: dimensions of matrices after all default auto- 
matic reductions in CPLEX 11, root relaxation time using CPLEX 11 Dual Simplex LP solver, and the remaining "elapsed time", after 
which CPLEX 11 is still at the root node of the monolithic formulation. For dives, the measurements are for the dive executed first in 
the given run. 



Instance 


Reduced LP (Monolithic) 


Reduced LP (Surface) 


Reduced LP (PeriodFixed) 


Reduced LP (DayFixed) 


compOl 


8593 x 6696 


29.01 s 


3807 x 2341 


2.01 s 


1484 x 1128 


0.14 s 


3984 x 4017 


5.33 s 


comp02 


34438 x 29559 


1058.47 s 


8629 x 5364 


9.11 s 


6605 x 5840 


2.80 s 


8089 x 10516 


16.76 s 


comp03 


30985 x 27240 


493.43 s 


7838 x 4966 


5.22 s 


5825 x 5136 


1.92 s 


7516 x 10230 


19.00 s 


comp04 


37570 x 33195 


590.15 s 


8097 x 5106 


1.51 s 


7367 x 6570 


2.29 s 


8991 x 14019 


19.53 s 


comp05 


19795 x 16652 


109.08 s 


9369 x 6946 


9.16 s 


2282 x 1827 


0.04 s 


5526 x 5474 


2.85 s 


comp06 


49360 x 43485 


1729.20 s 


10689 x 6645 


16.63 s 


9361 x 8442 


9.26 s 


13096 x 19977 


78.90 s 


comp07 


67901 x 60001 


6489.69 s 


13491 x 8108 


45.05 s 


12365 x 11300 


13.72 s 


16974 x 24940 


157.06 s 


comp08 


39935 x 35224 


899.81 s 


8643 x 5434 


2.63 s 


8240 x 7380 


6.02 s 


11316 x 17657 


42.66 s 


comp09 


36382 x 31951 


653.75 s 


8290 x 5333 


2.19 s 


7195 x 6390 


3.75 s 


8717 x 11881 


17.58 s 


complO 


51839 x 45759 


2087.70 s 


11244 x 6850 


21.58 s 


9628 x 8712 


8.97 s 


12388 x 20637 


83.90 s 


compll 


11747 x 8448 


28.86 s 


5699 x 3334 


2.93 s 


1367 x 960 


0.06 s 


5352 x 5896 


6.56 s 


compl2 


32893 x 27504 


480.57 s 


12570 x 8724 


27.25 s 


3990 x 3355 


0.30 s 


8114 x 8342 


10.84 s 


compl3 


39610 x 35188 


771.64 s 


8320 x 5329 


2.33 s 


8275 x 7410 


4.33 s 


10599 x 16795 


33.44 s 


compl4 


37785 x 32828 


866.63 s 


8685 x 5333 


8.10 s 


6905 x 6120 


2.23 s 


8151 x 12674 


13.19 s 



Fig. 8. Our results compared against the results of Lach and Liibbecke (2008a) and 
Miiller (2008). In all cases, CPLEX 11 is used with default parameters, no custom 
cuts, and CPU time normalised using the benchmark of Di Gaspero et al. (2007). 
Within 1 CPU unit, we use Surface and PeriodFixed. Within 10 or 40 CPU units, 
we use Surface2 with both PeriodFixed and DayFixed dives. 
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Compared to: Results of Liibbecke and Lach (using CPLEX 11) 
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7 Computational Experience 



The proposed heuristic has been implemented using ILOG Concert libraries 
in C++ (ILOG, 2007). ILOG CPLEX 11 is used both at Surface and in all 
dives. Some tests have also been carried out using ZIB SCIP, the present- 
best freely available integer programming solver developed at Konrad-Zuse- 
Zentrum fur Informationstechnik in Berlin (Achterberg, 2004); their results 
are available on request from the authors. The implementation has been 
evaluated on a single processor of a desktop PC equipped with two Intel 
Pentium 4 processors clocked at 3.20 GHz and with 4 GB of RAM, running 
Linux. To allow easier comparison with the results of Lach and Liibbecke 
(2008a) and entrants of the International Timetabling Competition 2007, 
time measurements have been normalised using the benchmark 3 by Di Gaspero 
et al. (2007). One CPU unit in Figure 8 corresponds to 780 seconds. Within 
the limit of 1 CPU unit, we spend 600 seconds solving Surface and then 180 
seconds in a single PeriodFixed dive. Within the limit of 10 CPU units, we 
spend 3420 seconds solving Surface2, then 180 seconds in a single Period- 
Fixed dive, and then 3600 seconds in a single DayFixed dive. Within the 
limit of 40 CPU units, we spend 3 hours solving Surf ace2, limited time in a 
number of PeriodFixed dives, and 5 hours in a single DayFixed dive. Full 
source code, logs, and solution files are available from the contact author 4 . 

Overall, the solver can produce feasible solutions and corresponding lower 
bounds for all instances from the International Timetabling Competition 
2007 within one or two CPU units. (It is only instances comp07 and compl2 
that take two CPU units.) Indeed, if some terms of the objective are ignored 
at Surface, then provided the relative weights of the other terms are not 
changed, the lower bound of the objective-restricted problem is also a valid 
lower-bound for the complete minimisation problem. However, within one 
CPU unit, alloted to solvers in the International Timetabling Competition 
2007, the progress is limited, mostly due to the sheer size of the LP in- 
stances involved, as outlined in Table 4. The gap 5 often remains as large 
as 95 %, if any solutions are found at all. The solver nevertheless contin- 
ues to make good progress, as suggested in Figure 8, and within ten and 
40 CPU units produces results with gaps comparable to those of Lach and 
Liibbecke (2008a). Compared to the results of Miiller (2008), the winner of 
Track 3 in the International Timetabling Competition 2007, the quality of 
the solutions obtained is only mediocre, no matter whether we consider the 
choice of the best result in ten separate runs of one CPU unit each equiva- 



3 http : / /www . cs . qub .ac.uk/ it c2007 /benchmarking/ 

4 http : //cs . nott . ac . uk/~jxm/ (Jakub Marecek) 

5 The gap is calculated using the formula of ILOG, i.e. 100(1 — LB/UB), where 
LB is lower bound and UB is cost of the best integer feasible solution. 
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lent to ten CPU units, as suggested by Lach and Lubbecke (2008a), or not. 
Except for the instance comp04, comp06, and comp08, we have not been 
able to improve upper bounds obtained by Miiller (2008) within 10 runs 
of 1 CPU unit within a single 40 CPU unit run. On comp04, however, we 
have closed the absolute gap to 1 (2.8%), a feat which seems rather diffi- 
cult to achieve using monolithic formulations. The possibility to perform 
a heuristic search whilst obtaining some information about its quality, per- 
haps together with the ease of implementation, hence represents two of the 
few, and perhaps the only two significant advantages of the proposed ap- 
proach over standard local search approaches. Nevertheless, it seems that, 
given this performance, such an approach could be well usable in practice, 
for example, in tightly constrained instances, where solutions delivered by 
local search heuristics are not good enough to accept, without proofs of 
their near-optimality. 



7.1 Remarks on the Instances of Linear Programming 



In this connection, it is also interesting to notice the properties of the lin- 
ear programming relaxations involved, summarised in Table 4. Generally, 
instances both at Surface and in all dives are considerably smaller than in- 
stances of Monolithic, and there seem to be a number of instances (comp06, 
comp07), where the linear programming solver based on barrier methods 
is, by orders of magnitude, faster. On difficult instances, such as comp06, 
the run time of the linear programming solver for Surf ace2 can be as much 
as four times as large as for Surface, and similarly many iterations can be 
performed per second later on. We do not have a good explanation for this 
behaviour. Compared to either Surface or Surface2, instances encountered 
in dives are very different. They are structurally more similar to weighted 
matching problems or multiple knapsack problems, often used to test in- 
teger programming solvers. PeriodFixed dives can often be solved to op- 
timality within a couple of minutes on instances from the International 
Timetabling Competition 2007. DayFixed dives often produce considerably 
better solutions within ten minutes on easier instances, but often do not pro- 
vide any feasible solutions within the first 30 minutes of the search on more 
difficult instances. Improving the numerical behaviour of linear program- 
ming solvers on the instances could perhaps lead to considerable improve- 
ments in performance. 
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8 Conclusions 



We have proposed a hybrid metaheuristic with an application in univer- 
sity course timetabling with multiple soft constraints. In addition to per- 
forming well on the particular problem of Udine Course Timetabling, this 
metaheuristic seems to provide a good starting point for the development 
of solvers for a number of similar problems with soft constraints with only 
minor modifications, including most timetabling problems with pattern pe- 
nalisation. The hybridisation we study might also be of more general inter- 
est. It is based on a mix of formulations of different sub-problems of the full 
problem, formulated in integer programming. The basic lesson learned is 
that the advantage of such hybridisation arises when the sub-problems are 
selected so as to have good properties that can help the integer program- 
ming solver. Specifically, the non-automatic reformulation of sub-problems 
is essential in order to gain full advantage. This is reminiscent of the dic- 
tum in local search, that the selection of neighbourhoods generally ought 
to be influenced by the issue of how effectively they can be implemented. 
Overall, this approach seems to present a promising direction in the devel- 
opment of hybrid metaheuristics. 

We believe that realistic problems often have many different terms in the 
objectives and that often these objective terms can be quite awkward to en- 
code in integer programming, requiring many extra variables. In this work 
on timetabling, we have seen an example of such 'messy' problems. We 
have seen that many different formulations, relaxations and abstractions 
can play a role. Furthermore, with control of their usage, results can be im- 
proved significantly with respect to direct use of the initial direct encod- 
ing. However, the number of possible ways to create relaxations and ab- 
stractions and link them together increases rapidly, and managing this can 
quickly become unyieldy. Hence, a long-term future goal is to provide in- 
telligent decision support to help the user keep track of the possibilties and 
semi-automatically discover good combinations. Our belief is that such a 
system would greatly help integer programming methods to achieve good 
results with less time and effort spent in 'tuning' the combinations of en- 
codings. 
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